systems which yield unique and exact solutions
Heliyon 4: e00691, 2018
https://doi.org/10.1016/j.heliyon.2018.e00691
Abstract
Many physical systems exhibit random or stochastic components which shape or even drive their dynamic behavior. The stochastic models and equations describing such systems are typically assessed numerically, with a few exceptions allowing for a mathematically more rigorous treatment in the framework of stochastic calculus. However, even if exact solutions can be obtained in special cases, some results remain ambiguous due to the analytical foundation on which this calculus rests. In this work, we set out to identify the conceptual problem which renders stochastic calculus ambiguous, and exemplify a discrete algebraic framework which, for all practical intents and purposes, not just yields unique and exact solutions, but might also be capable of providing solutions to a much wider class of stochastic models.
Supplementary Information
Original title and section headings
The presentation of this study somewhat follows a theme which I found to be of relevance for the subject matter being addressed. Specifically, the article's title along with the headings of its sections were chosen from key concepts marking the iconic motion picture "The Matrix" (1999) by The Wachowski Brothers. However, due to understandably strict formatting guidelines in the online journal Heliyon, I was kindly asked to replace these headings for the final published version. In the case an interested Reader would want to fully enjoy the intended presentational frame of this study, the following table lists the original headings.
|
A toy model not subject to the Itô-Stratanovich dilemma
Let us consider the following system of two first-order stochastic differential equations in x(t) and y(t): ˙x(t)=a1x(t)+a2x(t)y(t)+a3y(t),˙y(t)=b1y(t)+b2f(t), where all ai,bi∈Q:ai,bi≠0 are constants, and f(t) denotes an arbitrary stochastic processes, e.g. Gaussian white noise. The differential equation for y(t) encompasses the well-known stochastic model for Brownian motion, whereas x(t) describes a process multiplicatively driven by the latter. Although system (1) is not subject to the Itô-Stratanovich dilemma dealt with in the article, we will present here its treatment within the context of the proposed finite algebraic approach, as model (1) has far-reaching applications, for instance as effective stochastic model of neurons driven by a single multiplicative, i.e. conductance, synaptic noise source (e.g., see [1]).
Let us further consider f(t) as being a smooth function of t, then an explicit, formal solution of the above system (1) is given by x(t)=x(0)eI1(t)+a3t∫0dsI2(s)eI1(t)−I1(s),y(t)=y(0)eb1t+b2t∫0dsf(s)eb1(t−s), where I1(t)=t∫0ds(a1+a2I2(s)),I2(t)=y(0)eb1t+b2t∫0dsf(s)eb1(t−s). As an example, Fig. 1a illustrates for f(t)=sin(t) the numerical evaluation of this solution, and compares the latter to the numerical integration of the original system of differential equations (1).![]() |
Figure 1: Representative solutions of (1).
Compared are the numerical solutions of the original system (grey solid) and the explicit analytical solution (black solid and dotted). (a): Eqs. (2a) and (2b) for f(t)=sin(t); (b): Eqs. (4) and (9) for f(t)=const. Model parameters: a1=−1,a2=−1,a3=−1, b1=−1,b2=1, x(0)≡x0=0,y(0)≡y0=0, (a): f(t)=sin(t),(b): f(t)=1 |
The explicit analytic solution
Although Eqs. (2a) and (2b) provide, under the weak assumption of smooth and integrable f(t) in the interval [0,t], an explicit solution of system (1), a closed analytic form can, in general, not be obtained. A case which does, however, allow for such a solution is given if we assume f(t) as being piecewise constant. Specifically, let f(t)=ft be constant in the left-open finite interval (t,t+Δt], Δt>0. For notational simplicity, but without loss of generality, we will first consider the interval (0,Δt]. In this case, Eq. (2b) directly yields y(t)=y(0)eb1t+b2b1(eb1t−1)f0 ∀t∈(0,Δt].
The integration in Eq. (2a) is less straightforward and requires some more careful elaborations. Performing first the integration in (3a) and (3b), we obtain I1(t)=c1(eb1t−1)+c2t,I2(t)=c3(eb1t−1)+y(0), where c1=a2b2b21f0+a2b1y(0),c2=a1−a2b2b1f0,c3=y(0)+b2b1f0. With this, x(t) in Eq. (2a) takes the form x(t)=exp[c1eb1t+c2t]×{x(0)e−c1+a3c3b1Γ[−c2b1+1,c1;1,eb1t]+a3b1(y(0)−c3)Γ[−c2b1,c1;1,eb1t]}, where Γ[a,b;t0,t1] denotes the generalized double incomplete gamma function Γ[a,b;t0,t1]:=t1∫t0dtta−1e−bt.
The result (5) can be further simplified by utilizing various properties of the gamma function. To that end, we perform the coordinate transformation t→bt in (6), which yields for b≠0 Γ[a,b;t0,t1]=b−abt1∫bt0dtta−1e−t=b−a(γ[a,bt1]−γ[a,bt0]), where γ[a,t]:=∫t0dtta−1e−t=ta∞∑k=0(−1)kk!tka+k denotes the lower incomplete gamma function. Finally, using the recurrence relation γ[a,t]=(a−1)γ[a−1,t]−ta−1e−t, we can express (5) as x(t)=−a3a2+(x(0)+a3a2)exp[a2b21(b2f0+b1y(0))(eb1t−1)+1b1(a1b1−a2b2f0)t]−a1a3a2b1(a2b21(b2f0+b1y(0)))1b21(a1b1−a2b2f0)(γ[−c2b1,c1eb1t]−γ[−c2b1,c1])×exp[a2b21(b2f0+b1y(0))eb1t+1b1(a1b1−a2b2f0)t]. A representative example of this solution is illustrated in Fig. 1b.
Before proceeding, it is important to note that the original integral definition of γ[a,t] in Eq. (7) is only valid for a∈C with strictly positive Real[a]. However, its power series expansion can be shown to locally converge for all t∈C and a∈C/Z−:a≠0, where Z− denotes the set of all negative integer numbers. With this, Eq. (8) provides a valid recursion in the whole complex plane. In what follows, we will restrict to t∈Q and a∈Q/Z−, as well as γ[a,t]∈Q defined by the sum in (7), with the upper summation limit replaced by a sufficiently large yet finite number H chosen according to precision requirements or the characteristics of the underlying physical system. Although this choice, along with parameters and variables being assumed within Q and not R, might appear as irrelevant subtlety, it is essential for remaining within the discrete algebraic, i.e. finite, framework which forms the basis of the approach presented in the article. Moreover, each numerical or computational approach will naturally truncate expressions such as those occurring on the right-hand side of Eq. (7). Finally, we remark that, although γ[a,t] has singularities for a∈Z−, it can easily be shown by careful evaluation of the singular terms in the power expansion (7), that the difference between the two gamma functions occurring in the solution (9) yields finite values. Taking these arguments together, the solution for x(t) presented above and in the following sections, thus, must be considered as valid and is finite ∀a∈Q.
The recursive algebraic solution
As pointed out above, the solution of (1) for constant f(t)=f0∈Q and boundary values x(0),y(0)∈Q obtained in the previous section, Eqs. (4) and (9), is valid ∀t∈(0,Δt] with Δt>0. With an appropriate choice of boundary values, this solution can be easily generalized to arbitrary t>0. To that end, we assume that f(t) is a piecewise constant function with constant step width Δt, defined as f(t′)=ft=const ∀t′∈(t,t+Δt], with ft∈Q. In this case, the solution given in Eqs. (4) and (9) generalizes to y(t+s)=y(t)eb1s+b2b1(eb1s−1)ft and x(t+s)=−a3a2+(x(t)+a3a2)exp[a2b21(b2ft+b1y(t))(eb1s−1)+1b1(a1b1−a2b2ft)s]−a1a3a2b1(a2b21(b2ft+b1y(t)))1b21(a1b1−a2b2ft)(γ[−c2b1,c1eb1s]−γ[−c2b1,c1])×exp[a2b21(b2ft+b1y(t))eb1s+1b1(a1b1−a2b2ft)s], respectively, where s∈(0,Δt] and c1=a2b2b21ft+a2b1y(t),c2=a1−a2b2b1ft. As ft can take arbitrary values, Eqs. (11) and (12) provide a valid piecewise exact solution of (1) in intervals of length Δt even if f(t) describes a discrete-time stochastic process. Figure 2a visualizes a representative example in the case values of ft at different times t are drawn from a normal distribution, hence describes a discrete-time stochastic process which, in the statistical limit, resembles a Gaussian white noise process.
![]() |
Figure 2: Representative solutions of (1).
Compared are the numerical solutions of the original system (grey solid) solutions with the explicit analytical solution. (a): Eqs. (11) and (12) for f(t)=ft=const for t∈(t,t+Δt]) and recursive algebraic solution; (b): Eqs. (13) and (14) for fn=const ∀n∈N. Model parameters: a1=−1,a2=−1,a3=−1, b1=−1,b2=1, x(0)≡x0=0,y(0)≡y0=0, ft and fn were chosen from a normal distribution with mean 0 and standard deviation 0.4 |
Before proceeding towards a recursive algebraic form of the above solution, it is important to note that Δt used here plays a role which fundamentally differs from the notion of "step size" utilized in numerical integration schemes. As pointed out in the article, the latter serves as a mere quantitative measure of numerical discretization at the level of differential equations, and as such is subject to constraints to ensure not only the numerical stability of the solution, but also its numerical accuracy. However, in general, numerical integration can only deliver approximate solutions which typically accumulate numerical error. This is not the case in the solution presented above. Indeed, Eqs. (11) and (12) provide the exact analytic form of the solution of the system of differential equations (1) under the assumption (10), with Δt serving as a quantifier for the level of discretization of an analytically exact expression which is dictated solely by the "sampling rate" with which the input f(t) drives the system. Thus, in contrast to the notion of step size in numerical integration schemes, Δt is endowed with a direct link to physical reality. Recalling the fact that each physical observable is subject to constraints regarding its precision, each experimental observation can only deliver a finite and discrete set of values for a probed observable. However, such a set can always be expressed in a form similar to (10). In this sense, the solution in Eqs. (11) and (12) and, more generally, the approach presented here are consistent with the principles highlighted in Section 3 of the article, with Δt quantifying the limitations imposed by the process of measurement.
This independence of the numerical precision of our solution from Δt allows us to distance ourselves from the notion of "analysis" and, ultimately, the ideal yet unphysical construct of "continuity", and opens the possibility to deduce a recursive algebraic form of (11) and (12). Both of these equations are already recursive in the sense that the solution at t+s with s∈(0,Δt] depends solely on the solution at time t, and the constant value of f(t)=ft in the interval (t,t+Δt]. The values of the state variables x(s) and y(s) within open intervals s∈(t,t+Δt) thus become irrelevant due to the considered discretization, and can be discarded. With this, t acts as a mere label, and, without loss of generality, can be replaced by n∈N:n≥0. Equation (11) then takes the discrete, i.e. recursive, algebraic form yn+1=yneb1h+b2b1(eb1h−1)fn, where h:=Δt for notational convenience. Similarly, x(t) yields xn+1=−a3a2+(xn+a3a2)exp[a2b21(b2fn+b1yn)(eb1h−1)+1b1(a1b1−a2b2fn)h]−a1a3a2˜γ[b1,a2b2b21fn+a2b1yn,a1−a2b2b1fn]×exp[a2b21(b2fn+b1yn)eb1h+1b1(a1b1−a2b2fn)h], where ˜γ[a,b,c]=1abca(γ[−ca,beah]−γ[−ca,b])=H∑k=0(−1)kk!bkak−c(e(ak−c)h−1) with a,b,c∈Q, and h∈Q, H∈N denoting sufficiently small and large numbers, respectively, confined by experimental precision or other constraints imposed by physical reality. A representative example of the recursive algebraic solution is presented in Fig. 2b.
It is important to note that, due to their recursive nature, Eqs. (13) and (14) describe the exact incremental evolution of the system in finite and strictly positive steps h∈Q. In this sense, the obtained recursive algebraic solution is conceptually equivalent to a formulation in terms of differential equations, which describe the evolution of a given system in infinitesimal steps dt. Moreover, the algebraic nature of (13) and (14) ensures that the solution is not only unique, but also finite for arbitrary h and n, irrespective of the convergence properties for n→∞. With this, we next will recover the explicit algebraic form for xn and yn, thus move towards an exact solution of (1).
The explicit algebraic solution
Although Eqs. (13) and (14) appear to have a rather delicate mathematical form, both belong to the class of linear inhomogeneous recursions for which a whole host of techniques is readily available (e.g., see [2]). In order to apply the latter, we first introduce for notational convenience A:=eb1h,Bn:=b2b1(eb1h−1)fn and An:=exp[a2b21(b2fn+b1yn)(eb1h−1)+1b1(a1b1−a2b2fn)h],Bn:=−a3a2+a3a2exp[a2b21(b2fn+b1yn)eb1h+1b1(a1b1−a2b2fn)h]×{exp[−a2b21(b2fn+b1yn)]−a1˜γ[b1,a2b2b21fn+a2b1yn,a1−a2b2b1fn]}. With this, Eqs. (13) and (14) take the simpler forms xn+1=Anxn+Bn and yn+1=Ayn+Bn, respectively, with n≥0. In what follows, we will solve (18) and (19) utilizing what could be termed "operator approach", as both A and An act as operators on the respective state variables yn and xn at step n, thus evolving the system to step n+1.
First, we solve (19) for the state variable yn. It can easily be shown that repeated application of A yields yn=Any0+n−1∑i=0An−1−iBi for n∈N:n≥1, which, after inserting A and Bn defined in (16), takes the form yn=enb1hy0+b2b1(eb1h−1)˜Fn−1, where ˜Fn:=n∑i=0e(n−i)b1hfi. Equation (21) can be viewed as the explicit algebraic solution of the original differential equation for y(t) in (1) for piecewise constant inputs ft within a discrete finite mathematical framework, delivering the exact values of the state variable y(t) at discrete times t=nh. Interesting here is the occurrence of the term ˜Fn, which describes for each n the weighted sum over past inputs fi,0≤i<n. The presence of this term, however, is not surprising, as the original differential equation for y(t) describes a "filtered" version of the input f(t), i.e. a process which interlinks input values at different times through an exponentially decay with time constant b1. Although this, and the definition (22) of ˜Fn, suggest non-locality in time, ˜Fn is subject to the linear and, in general, inhomogeneous recursion ˜F0=f0,˜Fn+1=eb1h˜Fn+fn+1,n≥0. Thus, for arbitrary n, the value of ˜Fn only depends on its value ˜Fn−1 at the previous time step n−1 and the actual input fn at step n, so remains local in time.
In a similar fashion, successive application of Ai for i∈[0,n−1] in Eq. (18) yields the explicit expression xn=(n−1∏i=0Ai)x0+n−1∑i=1(n−1∏j=iAj)Bi−1+Bn−1 for the state variable xn with n∈N:n≥1. In contrast to yn, however, the evaluation of this expression is slightly less straightforward. After lengthy calculations (see section "Appendix: Exact algebraic solution of x(t)" below), we obtain an exact explicit algebraic solution for x(t) at discrete times t=nh,n∈N:n≥1 in form of xn=−a3a2+Mn−1ena1h(x0+a3a2−a1a3a2Pn−1), where Mn,n≥0 is given by the explicit expression Mn=exp[a2b2b21(eb1h−1)˜Fn−a2b2b1hFn+a2b1(e(n+1)b1h−1)y0], and the coefficients Pn by the inhomogeneous linear recursion P0=exp[a2b2b21f0+a2b1y0]˜γ[b1,a2b1b21f0+a2b1y0,a1−a2b2b1f0],Pn+1=Pn+exp[a2b2b1hFn+1+a2b2b21(1−b1h)fn+1+a2b1y0−(n+1)a1h]טγ[b1,a2b2b21(1−e−b1h)˜Fn+1+a2b2b21e−b1hfn+1+a2b1e(n+1)b1hy0,a1−a2b2b1fn+1], valid for all n≥0. Here, Fn denotes the linear (unweighted) sum over all previous inputs, i.e. Fn:=n∑i=0fi, obeying the linear recursive relation F0=f0,Fn+1=Fn+fn+1,n≥0.
The exact algebraic solution of (1) presented above in Eqs. (21) and (25) is interesting in various respects. Firstly, and most importantly, we note that, although both expressions are explicit in the initial state variables x0 and y0, the coefficients still contain recursive terms depending on the full history of the inputs fn in form of the auxiliary variables Fn and ˜Fn (see Eqs. (28) and (22), respectively). The presence of these terms is neither surprising nor conceptually at odds with results obtained when tackling the original system of differential equations (1) within the confines of differential calculus, i.e. an analytic approach. Integrating a differential equation, in fact the very concept of an "integral" in standard analysis, is by definition equivalent to an infinitesimally-paced summation over all functional values of the integrand. Although in the case of the differential system (1) with stochastic term f(t) such an integration is, in general, not possible in a mathematically rigorous sense without ambiguities, a finitely-paced summation within a discrete framework can be performed, as demonstrated above, leading to the presence of highly nonlinear yet finite and well-defined terms Mn and Pn in (26) and (27) containing the "integrated" history of the system.
Secondly, and somewhat surprisingly, the coefficients Mn and Pn allow, in the case of Mn, for a simple explicit and, in the case of Pn, for a linear recursive expression. Moreover, in both cases, for each given time step n, the defining expressions depend only on the current input fn and their linear and weighted sums Fn and ˜Fn, respectively, which by themselves follow simple recursive rules (Eqs. (28) and (22), respectively). Thus, values for xn and yn only depend on the current input and the value of a number of auxiliary variables at step n, i.e. are strictly local in time.
![]() |
Figure 3: Representative explicit algebraic solution of (1).
Compared are the numerical solutions of the original system (solid and dashed) solutions with the explicit algebraic solution (Eqs. (21) and (25); dots and triangles). (a): Constant input fn=const ∀n≥0 to systems with different intrinsic time constants with h=5; (b): fn=sin[(n+1)h]. Model parameters: (a): x0=1,y0=0,fn=0.5 ∀n≥0; slow dynamics (solid): a1=−0.01, a2=−1,a3=−1,b1=−0.1,b2=2; fast dynamics (dashed): a1=−10, a2=−1000,a3=−1000,b1=−1,b2=1; (b): a1=−1,a2=−1,a3=−1, b1=−1,b2=1,x0=0,y0=0 |
Finally, and as already detailed above, the parameter h has an interpretation which differs from the notion of "step size" in numerical solutions of differential equations. The latter is crucial for the accuracy of the numerical solution, especially when considering systems whose intrinsic dynamics is fast, in which case the integration step must be chosen small enough to faithfully capture dynamical properties. In contrast, the parameter h in Eqs. (21) and (25) is constraint by the sampling rate of the inputs fn, thus linked to experimental constrains or the resolution of the discrete-time stochastic process used in numerical approaches. To illustrate this point, Fig. 3a compares the case of constant input fn=const ∀n≥0 in two systems with fast and slow internal dynamics. Whereas in the former case, the numerical solution of the original differential equations requires a significantly smaller time step in order to capture the intrinsic dynamics (Fig. 3a, dashed), (21) and (25) yield precise solutions even if h is chosen much larger than the smallest time constant occurring in the original system. However, a more thorough investigation of this aspect of the presented discrete algebraic framework is required for a justifiable generalization, and lies outside the scope of this study.
Appendix: Exact algebraic solution of x(t)
For notational convenience, we define d1:=a2b2b21(eb1h−1−b1h),d2:=a2b1(eb1h−1),d3:=a1h,d4:=a3a2,d5:=a1a3a2b1,d6:=a2b2b21,d7:=a2b1,d8:=a1b1, as well as Ln:=d1fn+d2yn,Mn:=exp[n∑i=0Li],Nn:=d6fn+d7yn. Using (21), it can be shown that Ln allows for an explicit expression in terms of ˜Fn defined in (22), Ln=d1(˜Fn−˜Fn−1)+b2d2h˜Fn−1+d2y0enb1h for n≥1, and obeys the linear inhomogeneous recursion L0=d1f0+d2y0,Ln+1=eb1hLn+d1(fn+1−fn)+b2d2hfn,n≥0. In a similar fashion, the explicit form of Mn is given by Eq. (26), obeying the simple recursive relation M0=eL0,Mn+1=MneLn+1,n≥0. Finally, Nn is explicitly given by Nn=d6(˜Fn−˜Fn−1)+d7y0enb1h for n≥1, and follows the linear inhomogeneous recursion N0=d6f0+d7y0,Nn+1=eb1hNn+d6(fn+1−fn),n≥0. With this, An and Bn defined in (17) take ∀n≥0 the form An=eLn+d3,Bn=−d4+d4eLn+d3−b1d5eLn+Nn+d3˜γn, where ˜γn:=˜γ(a,b,c) with a=b1,b=d6(1−e−b1h)˜Fn+d6e−b1hfn+d7enb1hy0,c=a1−b1d6fn.
We first calculate xn for n=1. From Eq. (18) with n=0, or (24) with n=1, we immediately have x1=A0x0+B0=M0ed3x0−d4+d4M0ed3−b1d5M0eN0+d3˜γ0, which, with P0 defined in (27), is identical to (25) for n=1.
For n>1, we simplify the first term in (24): n−1∏i=0Ai=exp[n−1∑i=0Li+n−1∑i=0d3]≡Mn−1end3. Observing that n−1∏j=iAj=n−1∏j=0Aji−1∏j=0Aj≡Mn−1Mi−1e(n−i)d3, the second term in (24) takes the form (n−1∏j=iAj)Bi−1=−d4Mn−1Mi−1e(n−i)d3+d4Mn−1Mi−1eLi−1+(n−i+1)d3−b1d5Mn−1Mi−1eLi−1+Ni−1+(n−i+1)d3˜γi−1. Thus, for n>1, xn=−d4+Mn−1end3x0+d4Mn−1end3−b1d5Mn−1end3n−1∑i=0eLi+Ni−id3Mi˜γi, where we used the identity ed3n−1∑i=1eLi−1−id3Mi−1−n−1∑i=1e−id3Mi−1=1−e−(n−1)d3Mn−2, which is valid ∀n>1.
What remains is to simplify the sum occurring in the last term in xn above. The presence of ˜γn and its dependence on an arbitrary sequence fn renders the search for a closed-form or explicit expression difficult or perhaps impossible. However, observing that eLn+NnMn=exp[b1d6hFn+d6fn(1−b1h)+d7y0] and introducing ˜Pn:=exp[b1d6hFn+d6fn(1−b1h)+d7y0−nd3]˜γn ∀n≥0, we obtain xn=−d4+Mn−1end3x0+d4Mn−1end3−b1d5Mn−1end3n−1∑i=0˜Pi. Finally, defining Pn:=n−1∑i=0˜Pi and noting that the so-defined Pn obeys the recursive relation P0=˜P0,Pn+1=Pn+˜Pn+1,n≥0, we arrive at (25).
Dependency of the results on H
A representative example of the dependency of the algebraic solution on H is presented in Fig. 4, using a sinusoidal input as illustration. The parameters used in this example are identical to the system presented here (see Fig. 2b). For the purpose of simplicity, we restrict here to a consideration of H only for ˜γ[a,b,c], see (15). All other analytic functions, such as the exponential, are subject to the forced and exorbitant high precision goal of 10−100, applied throughout the study (in the article and the supplementary information presented here).
As can be seen, for H below 5, a significant error enters the algebraic solution for x(t). However, this is to be expected, as ˜γ[a,b,c] in these cases is approximated by the lowest orders in its power expansion. Specifically, we have ˜γ[a,b,c]=−ac(e−c−1)b−c/a for H=0, ˜γ[a,b,c]=ab−c/a(−b(ea−c−1)a−c−e−c−1c) for H=1, and ˜γ[a,b,c]=ab−c/a(b2(e2a−c−1)2(2a−c)−b(ea−c−1)a−c−e−c−1c) for H=2. Given the conceptual meaning of H as the upper limit of the power considered in the power expansion of analytic functions, it is then clear that such a heavily approximated solution must deviate from a solution which is obtained with the numerical precision goal of 10−100 used in the study.
![]() |
Figure 4: Representative examples of the solutions of the system of differential equations (1).
Compared are the numerical integration of the original system (solid) and the explicit analytical solution (25) for x(t) using a sinusoidal input fn=sin[(n+1)h] for different values of H. (a): H=0, (b): H=1, (c): H=2, (d): H=5, (e): H=10. The system's parameters are identical to those used for Fig. 2b above. |
Somewhat surprisingly, already a value of H=5 suffices to arrive at a solution which can be considered as indistinguishable from the "exact" (in the sense of a numerical evaluation with an exorbitantly high precision goal) solution. Despite this finding, in all numerical evaluations utilizing the algebraic solutions (recursive or explicit) throughout the study, a value of H=100 was used.
References[1] | Rudolph, M. and Destexhe A., Neural Comput. 15, 2577 (2003). |
[2] | Wilf, H. S., Generatingfunctionology (3rd edition, CRC Press, 2005). |